rm(list=ls())
library('dplyr')
library('tidyr')
library('VennDiagram')
library('grid')
library('futile.logger')
library("ggplot2")
library("gplots")
library('plotly')
library('tools')
library('cowplot')
library('UpSetR')
Import des 3 vcf obtenus après lancement du pipeline (vcf_12x, vcf_24x, vcf40x) et import du vcf de référence (vcf_ref).
vcf_12x = read.table("/home/elodie/Documents/Analysis/gold_12x_on_data_elodie_filtervarianttranches_2D_decomposed_normalized_uniq.vcf")
vcf_24x = read.table("/home/elodie/Documents/Analysis/gold_24x_on_data_elodie_filtervarianttranches_2D_decomposed_normalized_uniq.vcf")
vcf_40x = read.table("/home/elodie/Documents/Analysis/gold_40x_on_data_elodie_filtervarianttranches_2D_decomposed_normalized_uniq.vcf")
vcf_ref = read.table("/home/elodie/Documents/Elodie/Datas_ismael_ref/gold_on_data_elodie_ismael_decomposed_normalize_uniq.vcf")
plot1 = ggplot(nombre_variants, aes(fill=type, y=vcf, x=nombre)) +
geom_bar(position = "stack", stat = "identity") +
scale_fill_manual(values = c("#bebada", "#ffffb3", "#8dd3c7", "#170f0f")) +
labs(fill = "Type de variants", x = "Nombre de variants", y = "VCF") +
theme(panel.grid = element_blank(), plot.title = element_text(hjust = 0.5)) +
labs(title = "Nombre et type de variant par VCF")
ggplotly(plot1)
La distribution des types de variants semble équivalente pour chaque VCF. Cependant dans les 3 VCF générés par le pipeline on trouve des variants “Undetermined” mais pas dans le VCF de référence.
plot2 = ggplot(nombre_variants_par_gene, aes(fill=gene, y=vcf, x=nombre)) +
geom_bar(position="stack", stat="identity") +
scale_fill_manual(values = c("#fdc086", "#91bfdb", "#b2df8a")) +
labs(fill = "Nom des gènes", x = "Nombre de variants", y = "VCF") +
theme(panel.grid = element_blank(),plot.title = element_text(hjust = 0.5)) +
labs(title = "Nombre de variants par gène et par VCF")
ggplotly(plot2)
Le nombre de variants est identique dans chaque VCF pour les gènes UTS2 et LPL mais le nombre varie pour le plus grand gène EPHA5.
nombre_variants_par_gene$vcf_nom_genes = paste(nombre_variants_par_gene$vcf, nombre_variants_par_gene$gene, sep="_")
plot3 = ggplot(nombre_variants_par_gene, aes(fill=type, y=vcf_nom_genes, x=nombre)) +
geom_bar(position = position_fill(reverse=FALSE), stat="identity") +
scale_fill_manual(values = c("#bebada", "#ffffb3", "#8dd3c7", "#170f0f")) +
labs(title = "Nombre de variants par gène, par VCF et par type de variant") +
labs(fill = "Type de variants", x = "Nombre de variants", y = "VCF") +
theme(panel.grid = element_blank(), plot.title = element_text(hjust = 0.5))
ggplotly(plot3)
La répartition des types de variants par gène et par VCF est concordante. Aucune insertion dans UTS2 dans les 4 VCF.
plot4 = ggplot(nombre_type_snv[!is.na(nombre_type_snv$type), ], aes(fill = type, y = vcf, x = nombre))+
geom_bar(position = position_fill(reverse=FALSE), stat="identity") +
scale_fill_manual(values = c("#fc9272", "#9ebcda")) +
labs(title = "Type de SNV par VCF", x = "Type de substitutions", y = "Nombre de substitutions") +
labs(fill = "Type de SNV") +
theme(panel.grid = element_blank(), plot.title = element_text(hjust = 0.5))
plot4
En génétique humaine, il y a un environ 2 fois plus de transitions que de transversions. Sur le plot, on voit environ 2/3 de transitions et 1/3 de transversions ce qui est cohérent. Valeurs cohérentes dans tous les VCF.
plot5 = ggplot(nombre_type_snv[!is.na(nombre_type_snv$type_snv), ], aes(fill = type_snv, y = vcf_gene, x = nombre)) +
geom_bar(position = position_fill(reverse=FALSE), stat="identity") +
labs(title = "Répartition des types de SNV par VCF et par gène", x = "Nombre", y = "VCF et gène") +
scale_fill_manual(values = c("#a6cee3", "#1f78b4","#b2df8a", "#33a02c", "#fb9a99", "#e31a1c", "#fdbf6f", "#ff7f00", "#cab2d6", "#6a3d9a", "#ffff99", "#b15928")) +
labs(fill = "Type de SNV", x = "Nombre", y = "VCF") +
theme(plot.title = element_text(hjust = 0.5, size=14, face="bold"), plot.subtitle = element_text(hjust = 0.5)) +
theme(panel.grid = element_blank())
ggplotly(plot5)
Une répartition des types de SNV cohérente entre VCF pour chaque gène.
plot6 = ggplot(vcf_all %>% filter (type == "insertion" | type == "deletion"), aes(x = taille_delins, y = vcf)) +
geom_boxplot(aes(fill=type)) +
theme_classic() +
labs(x = "Taille des delins", y = "VCF", title = "Boxplots de la taille des insertions et délétions par VCF") +
theme(plot.title = element_text(hjust = 0.5, size=12, face="bold"), plot.subtitle = element_text(hjust = 0.5)) +
labs(fill = "Type de delins") +
scale_fill_manual(values = c("deletion" = "#bebada", "insertion" ="#ffffb3"))
plot6
La distribution des tailles des délétions est équivalente entre les 4 VCF. Les VCF 12x, 24x et 40x contiennent plus de valeurs aberrantes pour les insertions et délétions que le VCF de référence. La taille maximale des insertions du VCF de référence est plus petite que pour les 3 VCF obtenus avec le pipeline. L’interquartile est plus petit également pour la taille des insertions du VCF de référence ce qui signifie que les valeurs des insertions sont plus regroupées autour de la médiane. Il y a moins de disparité des valeurs pour le VCF de référence.
plot7 = ggplot(vcf_all %>% filter (type == "insertion"), aes(x = taille_delins, y = vcf)) +
geom_boxplot(aes(fill=gene)) +
theme_classic() +
labs(x = "Taille des insertions", y = "VCF", title = "Boxplots de la taille des insertions par VCF et par gène") +
theme(plot.title = element_text(hjust = 0.5, size=12, face="bold"), plot.subtitle = element_text(hjust = 0.5)) +
labs(fill = "Noms des gènes") +
scale_fill_manual(values = c("EPHA5" = "#fdc086", "LPL" ="#91bfdb", "UTS2" = "#b2df8a"))
plot8 = ggplot(vcf_all %>% filter (type == "deletion"), aes(x = taille_delins, y = vcf)) +
geom_boxplot(aes(fill=gene)) +
theme_classic() +
labs(x = "Taille des délétions", y = "VCF", title = "Boxplots de la taille des délétions par VCF et par gène") +
theme(plot.title = element_text(hjust = 0.5, size=12, face="bold"), plot.subtitle = element_text(hjust = 0.5)) +
labs(fill = "Noms des gènes") +
scale_fill_manual(values = c("EPHA5" = "#fdc086", "LPL" ="#91bfdb", "UTS2" = "#b2df8a"))
plot_grid(plot7, plot8)
plot9 = ggplot(vcf_12x_24x_40x, aes(x = as.numeric(vcf_12x_24x_40x$DP), y = vcf)) + geom_boxplot(aes(fill=gene)) + theme_classic() +
theme_classic() +
labs(x = "Profondeur de lecture des variants", y = "VCF", title = "Boxplots de la profondeur de lecture des variants par gène") +
theme(plot.title = element_text(hjust = 0.5, size=12, face="bold"), plot.subtitle = element_text(hjust = 0.5)) +
labs(fill = "Noms des gènes") +
scale_fill_manual(values = c("UTS2" = "#b2df8a", "LPL" ="#91bfdb", "EPHA5" = "#fdc086"))
plot9
## Warning: Use of `vcf_12x_24x_40x$DP` is discouraged.
## ℹ Use `DP` instead.
plot10 = ggplot(dp_12x, aes(x = DP, y = nombre)) +
geom_histogram(stat = "identity", fill = "#a6cee3") +
geom_vline(xintercept = mean(dp_12x$DP), color = "#fd8d3c", linetype = "dashed", size = 1) +
labs(x = "Profondeur de lecture", y = "Nombre de variants", title = "VCF 12x") +
theme(axis.text.x = element_text(hjust = 0.5), plot.subtitle = element_text(hjust = 0.5, size = 14, face = "bold")) +
theme(panel.grid = element_blank())
## Warning in geom_histogram(stat = "identity", fill = "#a6cee3"): Ignoring
## unknown parameters: `binwidth`, `bins`, and `pad`
## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
plot11 = ggplot(dp_24x, aes(x = DP, y = nombre)) +
geom_histogram(stat = "identity", fill = "#1d91c0") +
geom_vline(xintercept = mean(dp_24x$DP), color = "#fd8d3c", linetype = "dashed", size = 1) +
labs(x = "Profondeur de lecture", y = "Nombre de variants", title = "VCF 24x") +
theme(axis.text.x = element_text(hjust = 0.5), plot.subtitle = element_text(hjust = 0.5, size = 14, face = "bold")) +
theme(panel.grid = element_blank())
## Warning in geom_histogram(stat = "identity", fill = "#1d91c0"): Ignoring
## unknown parameters: `binwidth`, `bins`, and `pad`
plot12 = ggplot(dp_40x, aes(x = DP, y = nombre)) +
geom_histogram(stat = "identity", fill = "#253494") +
geom_vline(xintercept = mean(dp_40x$DP), color = "#fd8d3c", linetype = "dashed", size = 1) +
labs(x = "Profondeur de lecture", y = "Nombre de variants", title = "VCF 40x") +
theme(axis.text.x = element_text(hjust = 0.5), plot.subtitle = element_text(hjust = 0.5) ) +
theme(panel.grid = element_blank())
## Warning in geom_histogram(stat = "identity", fill = "#253494"): Ignoring
## unknown parameters: `binwidth`, `bins`, and `pad`
plot_grid(plot10, plot11, plot12, ncol = 3, nrow = 1)
On ne fera pas d’exploration de la VAF ni de la profondeur DP pour le VCF de référence car il n’y a pas d’explication aux valeurs de DP, AD, ADALL dans le VCF. Les valeurs ne sont pas cohérentes avec ce qui est observé dans le bam de référence associé au VCF dans IGV.
cat("VAF minimale dans VCF 12x :", min(vaf_12x$vaf), "%", "\n")
## VAF minimale dans VCF 12x : 16.7 %
cat("VAF minimale dans VCF 24x :", min(vaf_24x$vaf), "%", "\n")
## VAF minimale dans VCF 24x : 14.3 %
cat("VAF minimale dans VCF 40x :", min(vaf_40x$vaf), "%", "\n")
## VAF minimale dans VCF 40x : 16.3 %
plot13 = ggplot(vcf_vaf, aes(x = nombre_variants, y = cat_vaf , fill = vcf)) +
geom_bar(position = position_fill(reverse=FALSE), stat="identity") +
labs(title = "Répartition du nombre de variants par catégorie de VAF", x = "Nombre de variants", y = "Catégorie de VAF") +
scale_fill_manual(values = c("#a6cee3", "#1d91c0", "#253494")) +
labs(fill = "VCF") +
theme(plot.title = element_text(hjust = 0.5, size=14, face="bold"), plot.subtitle = element_text(hjust = 0.5)) +
theme(panel.grid = element_blank()) +
theme_classic()
ggplotly(plot13)
variants_list = list(vcf_12x = unique(vcf_12x$variants), vcf_24x = unique(vcf_24x$variants), vcf_40x = unique(vcf_40x$variants), vcf_ref = unique(vcf_ref$variants))
plot14 = upset(fromList(variants_list), order.by = "freq", text.scale = 2, point.size = 5, sets.bar.color = "#999999", main.bar.color = "#fc8d62")
plot14
plot15 = display_venn(variants_list, category.names = c("vcf_12x" , "vcf_24x" , "vcf_40x", "vcf_ref"),
fill = c("#a6cee3", "#1d91c0", "#253494", "#E69F00"), cat.cex = 1.2,
cat.fontface = "bold",
cat.default.pos = "outer",
cat.dist = c(0.1, 0.1, 0.1, 0.1))
plot16 = ggplot(cat_vcf_all_nombre, aes(fill=type, x=nombre, y=categorie)) +
geom_bar(position = "stack", stat="identity") +
scale_fill_manual(values = c("#bebada", "#ffffb3", "#8dd3c7", "#170f0f")) +
labs(title = "Nombre de variants par catégories et par type de variants") +
labs(fill = "Type de variants", x = "Nombre de variants", y = "Catégories de VCF") +
theme(panel.grid = element_blank(), plot.title = element_text(hjust = 0.5))
plot17 = ggplot(cat_vcf_all_genes, aes(fill=gene, x=nombre, y=categorie)) +
geom_bar(position = "stack", stat="identity") +
scale_fill_manual(values = c("UTS2" = "#b2df8a", "LPL" ="#91bfdb", "EPHA5" = "#fdc086")) +
labs(title = "Nombre de variants par catégories et par gènes") +
labs(fill = "Type de variants", x = "Nombre de variants", y = "Noms des gènes") +
theme(panel.grid = element_blank(), plot.title = element_text(hjust = 0.5))
plot_grid(plot16, plot17, ncol = 1, nrow = 2)
variants_list_12x_24x_40x = list(vcf_12x = unique(vcf_12x$variants), vcf_24x = unique(vcf_24x$variants), vcf_40x = unique(vcf_40x$variants))
plot18 = display_venn(variants_list_12x_24x_40x, category.names = c("vcf_12x" , "vcf_40x", "vcf_24x"),
fill = c("#a6cee3", "#1d91c0", "#253494"), cat.cex = 1.2,
cat.fontface = "bold",
cat.default.pos = "outer",
cat.dist = c(0.1, 0.1, 0.1))
recall (sensibilité), précision (spécificité) et score F1. Ce sont des métriques qui permettent de voir si un modèle est performant. Le score F1 combien recall et précision Plus le recall est élevé, plus le modèle repère des positifs. Plus la précision est élevée, moins le modèle se trompe sur les positifs (moins il y a de faux positifs). Plus le score F1 est élevé, plus le modèle est performant.
happy_12x_summury = read.table("/home/elodie/Documents/Analysis/Happy_12x/Happy_12x/ref-12x.summary.csv", header=TRUE, sep=',')
happy_24x_summury = read.table("/home/elodie/Documents/Analysis/Happy_24x/Happy_24x/ref-24x.summary.csv", header=TRUE, sep=',')
happy_40x_summury = read.table("/home/elodie/Documents/Analysis/Happy_40x/Happy_40x/ref-40x.summary.csv", header=TRUE, sep=',')
plot_roc_snp = plot_data(roc_snp, is.PR=TRUE)
## Warning: `aes_string()` was deprecated in ggplot2 3.0.0.
## ℹ Please use tidy evaluation idioms with `aes()`.
## ℹ See also `vignette("ggplot2-in-packages")` for more information.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
plot_roc_indel = plot_data(roc_indel, is.PR=TRUE)
plot_grid(plot_roc_snp, plot_roc_indel, nrow=2, ncol=1)